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Background 


Currently,  in  the  U.S.  there  is  a  movement  to  construct  one  NWP  model  (NWS, 
Navy,  Air  Force  -  other  partners  include  NASA  and  DOE).  This  National 
Board  (NUOPC=National  Unified  Operational  Prediction  Capability)  aims  to 
develop  a  new  model  that  is: 

1 .  Highly  scalable  on  current  and  future  computer  architectures 

2.  Global  model  that  is  valid  at  the  meso-scale  (i.e.,  non-hydrostatic) 

3.  Applicable  to  medium-range  NWP  and  decadal  time-scales 


Motivation 


Our  goal  is  to  construct  numerical  methods  for  non-hydrostatic  mesoscale  and 
global  atmospheric  models  (for  NWP  applications). 


To  verify  our  numerical  methods  we  have  built  a  modeling  framework  with  the 
following  capabilities: 

1 .  Highly  scalable  on  current  and  future  computer  architectures  (exascale 
computing:  this  means  CPUs  and  GPUs) 

2.  Flexibility  to  use  a  wide  range  of  grids  (e.g.,  statically  and  dynamically 
adaptive) 

3.  Model  that  is  accurate,  robust,  and  fully  conservative 


Talk  Summary 


1 .  New  models  need  to  exploit  available  computers 

2.  Numerical  methods  in  new  GFD  models 

3.  Desired  properties  of  our  future  GFD  models 

4.  NUMA  Model 


Talk  Summary 


1 .  New  models  need  to  exploit  available  computers 

•  Clock-speeds  no  longer  getting  faster;  vendors  are  just  giving 
us  more  CPUs 

•  From  Terascale  to  Petascale/Exascale  Computing 

•  10  of  Top  500  are  already  in  the  Petascale  range 

•  3  of  top  10  are  GPU-based  machines 

2.  Numerical  methods  in  new  GFD  models 

3.  Desired  properties  of  our  future  GFD  models 

4.  NUMA  Model 


Example  of  Linear  (Perfect)  Scalability 


x  (m) 


512  1024  2048  4096  8192  16384 

number  of  processors 


NUMA-CG  Simulation  with  16  Million  Grid  Points 


Talk  Summary 


New  models  need  to  exploit  available  computers 

Numerical  methods  in  new  GFD  models 

Time-Integration  is  important  (e.g.,  explict,  fully-implicit, 
semi-implicit) 

Spatial  Discretization  methods  is  how  we  are  able  to  take 
advantage  of  Parallel  computers  (i.e.,  domain  decomposition 
of  the  physical  grid) 

Desired  properties  of  our  future  GFD  models 
NUMA  Model 


Time  Integration 


Explicit  methods 
Fully-implicit  methods 
Implicit-Explicit  (IMEX)  methods 

-  Linear  Multi-step  Methods 

-  Multi-stage  Methods 

-  Multi-rate  Methods 


Multi-step/Multi-stage  IMEX  Methods 


Let’s  write  the  governing  equations 
as 


dq 

dt 


=  S(q) 


If  we  knew  the  linear  operator  L 
(containing  the  fastest  waves  in  the 
system),  then  we  could  write 


^  =  {S(q)  -  <5s;L(q)}  +  Ssl  [L(q)] 


Discretizing  by  a  Kth  order  time- 
integrator  yields 


> 


q„  =q  +  AL(q„) 


(/-AL)q„  =  q 


> 


Ax  =  b 


IMEX  Methods 

(Important  Properties  of  Schur  Complement) 
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The  Schur  system  is  clearly  smaller  than 
the  original  (No  Schur)  system  (NxN 
instead  of  3Nx3N  for  2D  SWE). 

Equally  important  is  that  the  Schur 
System  is  better  conditioned  than  the 
original  system.  This  means  that  fewer 
iterations  are  required  by  an  iterative 
solver  to  reach  convergence. 

Key  Point:  Eigenvalues  of  the  Schur 
form  lie  in  the  region  where  GMRES  will 
perform  well. 
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2D  Euler  Equations 


Multi-Rate  Methods 


3  6 

W-E  [1 000  Km] 


Adaptive  Mesh  Refinement 
(Space) 


Multi-rate  Time-Integration 
(Space-Time) 


Spatial  Discretization  Methods 


Element-based  Galerkin  Methods 

-  Continuous  Galerkin 

-  Discontinuous  Galerkin 


Element-based  Galerkin  (EBG)  Methods 

(Definition  and  Examples) 


An  element  is  chosen  to  be  the  basic  building-block  of  the  discretization 
and  then  a  polynomial  expansion  is  used  to  represent  the  solution  inside 
the  element 


Element-based  Galerkin  Methods 
(N=3  Basis  Function  Expansion) 


Spatial  Discretization 


Primitive  Equations: 


dq 

dt 


+  V-F  =  S(q) 


Approximate  the  solution  as: 

-  Interpolation  O(N) 

Write  Primitive  Equations  as: 


mn 

c1n  ~  2>,  tfi 


i= 1 


Mn 

dt 


+  V-F; 


N 


Weak  Problem  Statement:  Find 
qN  e  E(Q)  V  y/  e  S 


■  2  =  {^e//‘(Q) 
E  =  {i//eL2(Q) 


L  ^ 


f(?Af)  SN—S(qN) 

-SN=e 

■.¥eP,,(nmi)  (CG) 
W  e  ^Af(Oe)V£2eJ  (DG) 

)d£le  =  0 


such  that  (Integration  0(2N) ) 


Spatial  Discretization 


Integral  Form: 

J*^  \l/R(qN)d£l( 

Matrix  Form: 

(D!f) 

Communicator: 

C  (L(q)) 

Where  each  matrix  is:  M(e)—  f  ur  ur  WO 

V1  ij  jQ  r  iY  ja^^e 

DT=iiJy/.y/jdne 

Integration  0(2N) 


For  DG:  C(L,(9))  =  L,(9)  +  (m£  f  F<* 


K  = 


i  =  Jrn  VWjdT 


For  CG: 


C{Li(q))  =  S(G(q)) 


Integral  Form: 

Matrix  Form: 
Communicator: 


Spatial  Discretization 

L  YR(qN)d&e  =0 

dq '  ■  (d^  )'  Fw-  S(‘ >=  0 


LAq)  =  Mf^J 


dt 


c  (L(q)) 


Where  each  matrix  is: 


M 


t=  L  < 


D®)=  L V V.Vj^K 


For  DG:  C{L,(qj)  =  Li(q)  +  (M^j  f  F<* 

For  CG:  C(L,(q))  =  S(G(q)) 


Integration  0(2N) 


-»  M~  = 


£  =  j 


■»/ 


q,  =GW,e>)  (he) 


Spatial  Discretization 


Integral  Form: 

J*^  \l/R(qN)d£l( 

Matrix  Form: 

=  (DD 

Communicator: 

C  (L(q)) 

Where  each  matrix  is:  M(e)—  f  ur  ur  WO 

V1  ij  jQ  r  iY  ja^^e 

DT=iiJy/.y/jdne 

Integration  0(2N) 


For  DG:  C(L,(9))  =  L,(9)  +  (m£  f  F<* 


M!/  =  ln¥i¥jdT 


For  CG:  C(L,^))  =  S(Gfo)) 


q\e)  =s(<7/)  1 - K*» 


Domain  Decomposition:  Adjacency  Matrices 
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Geometry 
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DG 

Quadrilaterals 
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CG  Adjacency 


DG  Adjacency 


Talk  Summary 


New  models  need  to  exploit  available  computers 

Numerical  methods  in  new  GFD  models 

Desired  properties  of  our  future  GFD  models 

E.g.,  Conservation,  Scalability,  High-order  Accuracy, 
Adaptivity 

NUMA  Model 


Desired  Properties  of  Future  GFD  Models 


Conservation  -  Conservation  of  Mass  and  Energy  are  absolute  musts; 
what  else  should  we  conserve? 

Scalability  -New  models  must  be  highly  scalable  because  we  will 
continue  to  get  more  processors 

High-Order  Accuracy  -  Accuracy  is  important,  of  course,  but  how  do  we 
measure  this  and  what  order  accuracy  is  sufficient?  This  question  is 
coupled  to  the  accuracy  of  the  physics,  data  assimilation,  etc.  From  the 
standpoint  of  scalability,  high-order  is  good  (hp  methods  =  on-processor 
work  is  large  but  the  communication  footprint  is  small).  This  is  also  a 
good  strategy  for  exploiting  MPI/Open  MP  Hybrid  and  CPU/GPU 
paradigms. 

Adaptivity  -  Adaptive  methods  have  improved  tremendously  in  the  past 
decade  and  it  may  offer  an  opportunity  to  solve  problems  not  feasible  a 
decade  ago  but  we  need  to  identify  these  applications  (e.g.,  hurricanes, 
storm-surge  modeling,  clouds?). 
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NUMA-CG  and  NUMA-DG  with  16  Million  Grid  Points 
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Adaptivity 

(Muller,  Behrens,  Giraldo,  Wirth  2011) 
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Rising  Thermal  Bubble  with  Nonhydrostatic  DG  Model 


Talk  Summary 


New  models  need  to  exploit  available  computers 
Numerical  methods  in  new  GFD  models 
What  should  we  aim  for  in  our  new  models 
NUMA  Model 

Design  Philosophy 
Governing  Equations 
Grids 
Results 


Design  Philosophy 


Unified  Dynamics 

•  All  limited-area  models 
are  nonhydrostatic. 
Resolutions  of  global 
models  are  approaching 
the  nonhydrostatic  limit 
(~10  km). 

•  Both  limited-area  and 
global  models  utilize  the 
same  equations. 

•  Engineer  a  common 
dynamical  core 
(DyCore)  for  both 
models,  then  change 
grids  and  boundary 
conditions. 


•  Unified  Numerics 

•  CG  is  more  efficient  for 
smooth  problems  at 
low  processor  counts. 

•  DG  is  more  accurate 
for  problems  with 
sharp  gradients  and 
more  efficient  at  high 
processor  counts. 

•  Both  EBGs  utilize  a 
common  mathematical 
arsenal. 

•  NUMA  allows  the  user 
to  choose  either  CG  or 
DG  for  the  problem  at 
hand. 


Unified  Code 

Code  is  modular , 
with  a  common  set 
of  data  structures. 

New  time- 
integrators,  grids, 
basis  functions, 
physics,  etc.  may 
be  swapped  in  and 
out  with  ease. 

Code  is  portable: 
Successfully 
installed  on  Apple, 
Sun,  Cray,  and 
IBM. 


Governing  Equations 

(Unified  Global/Mesoscale  Equations) 

Beginning  with  the  Equations: 

^  +  V*(pu)  =  0 
ot 

—  +  u*  Vu  +— VP  +  gr  +  /(?xu)  =  0 
dt  p 

— +u* V0=O 
dt 


A 


pR0  ^ 


rA  ) 


We  introduce  the  General  Splitting: 


p(x,  y,  z,  t)  =  p0  (x,  y,z)  +  p  \x,  y,  z,  t) 
0(x,  y,  z,  t)  =  90  (x,  y,  z)  +  9  '(x,  y,z,t) 
P(x,y,z,t )  =  P0(x,y,z)  +  P'(x,y,z,t ) 


Governing  Equations 

(Unified  Global/Mesoscale  Equations) 

The  Reference  Fields  satisfy  certain  conditions,  e.g.,: 

f  •  VP0  =  -p0g 

— VP0+£r  +  /(rxu)  =  0 
Po 

The  following  General  Equations  satisfy  such  conditions: 

^  +  V.(pu)  =  0 
ot 

^-  +  u  •  Vu  + — (HVP0  +  VP'+  p'gr )  =  - /(rxu) 

^  +  u*V6»  =  0 
dt 


Mesoscale  Modeling  Mode 


Grids 

Global  Modeling  Mode 


Global  Modeling  Mode 
(Icosahedral) 


Telescoping  Grid 


ITCZ  Grid 


Results:  Mesoscale  Mode 

(3D  Linear  Hydrostatic  Isolated  Mountain) 


•  Flow  ofU=20  m/s  in  an  isothermal  atmosphere. 

•  LH  Mountain:  Solid  of  revolution  of  Witch  of  Agnesi:  Mountain  height  =  1  m 
with  radius  10  km. 

•  Absorbing  (sponge)  boundary  condition  implemented  on  lateral  and  top 
boundaries. 


Grid 


Geometry 


3D  Linear  Hydrostatic  Isolated  Mountain 

(T=l  hour) 


80  100  120  140  160 
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3D  Linear  Hydrostatic  Isolated  Mountain 

(T=l  hour) 


x  1CT3 


Normalized  L  Norm 


Results:  Mesoscale  Mode 

(Balance  Initial  State  and  Baroclinic  Instability) 


Balanced  Initial  State  Baroclinic  Instability 

(Hydrostatically  and 

Geostrophically 

balanced) 


Results:  Global  Mode 

(Planetary  Acoustic  Wave  Propagation) 


t  =  0.083  hours 


Theory=347  m/s, 

NUMA  model=347  m/s, 
NIC  AM  model  =  338  m/s 
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Results:  Global  Mode 

(Inertia-Gravity  Wave  Propagation  N=0.01,  T=48  hours) 


longitude  (degrees) 
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Theory=32  m/s, 
NUMAmodel=33  m/s, 
NIC  AM  model  =  33  m/s 


Rain-rate  (mm/h) 


NUMA2D  Simple  Moist  Physics 

(Mesoscale) 


qc  (g/kg) 


A  Multitude  of  Challenges  Remain 


Continuous  and  Discontinuous  Galerkin  methods  are  good  choices  for 
hydrostatic  and  non-hydrostatic  atmospheric  models. 

The  NUMA  dynamical  core  is  quite  mature: 

-  3D  and  MPI  (Kelly  and  Giraldo  JCP  201 1) . 

-  Can  use  either  CG  or  DG  methods. 

-  Contains  a  suite  of  IMEX  time-integrators  (Giraldo  et  al.  SISC  2011). 

-  Has  been  verified  on  a  variety  of  limited-area  and  global  tests. 

The  NUMA  physics  is  under  development: 

-  Simple  sub-grid  scale  parameterization  has  been  added  to  NUMA2D 
(Gabersek  et  al.  MWR  201 1). 

-  Variational  Multi-scale  (VMS)  method  is  being  added  to  NUMA2D  with 
physics  for  positivity  of  tracers  (Marras  et  al.  JCP  201 1). 

-  Simple  physics  being  added  to  NUMA3D  along  with  VMS  algorithm. 


A  Multitude  of  Challenges  Remain 


Future  Projects: 

-  Conduct  a  variety  of  verification  and  validation  simulations. 

-  Add  IMEX  methods  to  NUMA-DG  (currently  only  in  NUMA-CG  part). 

-  Study  and  Implement  Multi-Rate  (Patrick  Mugg  Thesis). 

-  Explore  new  Riemann  solvers  for  NUMA-DG  (Maria  Lukacova) 

-  Develop  adaptive  methods  for  unified  limited-area/global  modeling 
simulations  (work  in  progress  by  Kopera  and  Gopalakrishnan,  NPS) 

-  NUMA  is  being  coupled  with  PETSc  -  ANL  (DoE  ASCR  program)  has 
chosen  NUMA  as  its  flagship  application. 


